【汎用性抜群】尤度比検定をわかりやすく解説します

こんにちは、青の統計学です。

尤度比検定とは、汎用性の高い統計モデルの検定です。

その汎用性の高さは、サンプル数が十分大きい時には、尤度比検定統計量の対数に2をかけたものがカイ２乗分布に従う性質にあります。

行列を使わず解説したコンテンツは以下になります。

【python】尤度比検定で統計モデルの比較をしよう|統計的仮説検定

尤度比検定(likelihood ratio test)

まず、専門的な用語抜きに説明すると、尤度比検定とは二つのモデルのうち、観測データをよりよく説明するのはどちらだろうか、という問いに答えるための統計的な手法です。

ここで一方のモデルの尤もらしさをもう一方のモデルの尤もらしさで割ってあげることで、統計量を作ります。「尤もらしさ」が「最大尤度」になります。

尤度についておさらいしたい方はこちらをご覧ください。

【尤度って？】尤度関数と最尤推定量の解説と例題

繰り返しになりますが、尤度比検定（Likelihood ratio test, LRT）は統計学において、2つの統計モデルの適合度を比較するために使用される手法です。

この検定は、ネストされたモデル（すなわち、一方のモデルが他方のモデルに含まれるパラメータの制約版である場合）の適合度を比較する際によく使われます。

では基本的なアプローチを説明します。

モデルの設定

完全モデル（制約なしモデル）と縮小モデル（制約ありモデル）を設定します。

完全モデルはパラメータがより多く、一方の縮小モデルはそのパラメータの一部に制約（例えば、あるパラメータが0であるなど）が課されます。

完全モデル

$$\mathbf{y}=\mathbf{X}\beta+\epsilon$$

・$\mathbf{y}$は$n×1$の目的変数ベクトルです。

・$\mathbf{X}$は、$n×(p+1)$の設計行列です。ちなみに、説明変数の値を含んでします。1列目はすべて1でバイアス項$\beta_0$のためのものです、

・$\beta$はパラメータベクトルで、$(p＋1)×1$です。

・$\epsilon$は$n×1$の誤差項ベクトルです。

縮小モデル

もし $\beta_2=0$ という制約を設ける場合、このパラメータをモデルから除外します。行列 $\mathbf{X}$ から該当する列 $x_2$ を削除し、パラメータベクトル$\beta$からも$\beta_2$ を除きます。この場合のモデルは次のように表されます：

$$\mathbf{y}=\mathbf{X}_{-2}\beta_{-2}+\epsilon$$

・$\mathbf{X}_{-2}$は、$\mathbf{X}$から$\beta_2$に対応する列を除いた$n×p$の設計行列です。

・$\beta_{-2}$は$\beta_2$を除いた$p×1$のパラメータベクトルです。

尤度比検定統計量の作成

尤度比検定では、これら2つのモデルの尤度を比較します。

尤度比 $\Lambda$ は次のように定義されます

$$\Lambda=\frac{L(\beta_{-2}|\mathbf{X}_2,\mathbf{y})}{L(\beta|\mathbf{X},\mathbf{y})}$$

分子が縮小モデルの尤度で、分母が完全モデルの尤度ですね。

さて、尤度比から尤度比検定統計量を作りましょう。

$${\lambda=-2log\Lambda=-2log\frac{L(\beta_{-2}|\mathbf{X}_2,\mathbf{y})}{L(\beta|\mathbf{X},\mathbf{y})}}$$

尤度比検定量$\lambda$は、大標本（大きなサンプルサイズ）の下で、自由度が完全モデルと縮小モデル間のパラメータ数の差に等しいカイ二乗分布に従います。

この性質は大標本の理論の一部ウィルクスの定理に基づいています。

$${\lambda\sim \chi^2_{df}}$$

ここで、$df$は自由度で、完全モデルと縮小モデルの間のパラメータ数の差です。

カイ二乗分布は、F分布やt分布と同様に標本に関する分布と呼ばれています。

こちらについては、本記事の末尾に補足をしております。

CODE｜一般化線形モデル（GLM）における尤度比検定

ちょっとした検定を行います。

ポアソン回帰モデルを用いて、ある説明変数が応答変数に対して有意な影響を持つかどうかを尤度比検定で評価します。

ポアソン回帰モデルは、一般化線形モデルの一つです。

詳しくはこちらのコンテンツをご覧ください。

【GLM】一般化線形モデルを解説|ポアソン回帰

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import chi2

# データ生成
np.random.seed(0)
data = np.random.poisson(lam=30, size=1000)
covariate = np.random.normal(size=1000)
data_frame = pd.DataFrame({'Response': data, 'Covariate': covariate})

# 完全モデル（Covariateを含む）
full_model = smf.glm('Response ~ Covariate', data=data_frame, family=sm.families.Poisson()).fit()

# 縮小モデル（Covariateを除外）
reduced_model = smf.glm('Response ~ 1', data=data_frame, family=sm.families.Poisson()).fit()

# 尤度比検定
lr_stat = reduced_model.deviance - full_model.deviance
lr_df = full_model.df_model - reduced_model.df_model
lr_pvalue = chi2.sf(lr_stat, lr_df)  # カイ二乗分布からのp値を計算

print(f"LR statistic: {lr_stat}")
print(f"p-value: {lr_pvalue}")
print(f"Degrees of freedom: {lr_df}")

LR statistic: 1.9353660252947975

p-value: 0.16417256842907785

Degrees of freedom: 1

次元を一つ落としただけなので、自由度は1です。

p値の部分を補足をすると、カイ二乗分布の生存関数（1 – CDF）を計算しています。

これは、検定統計量がカイ二乗分布で観察された値よりも大きくなる確率を返し、これによってp値が得られます。